MatSqrtInverse <- function(A) {
    ##  Compute the inverse square root of a matrix
    ei <- eigen(A)
    d <- pmax(ei$values, 0)
    d2 <- 1/sqrt(d)
    d2[d == 0] <- 0
    ## diag(d2) is d2 x d2 identity if d2 is scalar, instead we want 1x1 matrix
    ei$vectors %*% (if (length(d2)==1) d2 else diag(d2)) %*% t(ei$vectors)
}

## Convert list to data frame
ConvertSETypes <- function(s){
    Vhat <- names(s)
    do.call(rbind,
            lapply(Vhat, function(v) data.frame(Vhat=v, dof=s[[v]],
                                                stringsAsFactors=FALSE)))
}

FindMinInterval <- function(f, x, tol=0.01, maxlength=30) {
    ## Construct confidence set by inverting a (non-randomized) test "f". f
    ## rejects H0:theta=x0 if f(theta0)=TRUE, and accepts H0 if f(theta0)=FALSE.
    ## User supplies initial interval "x". Numerical tolerance given by "tol".
    ## Method assumes that confidence set is an interval. In case interval ends
    ## up being longer than maxlength, return the interval c(0, maxlength)

    FinerGrid <- function(x){
        ## make finer grid of xs
        x1 <- rep(x, each=2)
        (x1[-1] + x1[-length(x1)]) / 2
    }

    ## Produce vector x of length 4 such that f(x)=c(1,0,0,1)
    FindZero <- function(f, x){
        ## in case interval is too long
        if (x[length(x)]-x[1] > maxlength) return(maxlength)

        fx <- sapply(x, f)
        gs <- x[length(x)]-x[1]         # grid size

        if (fx[1]==0) {                 # make left grid endpoint smaller
            FindZero(f, c(x[1]-gs/2, x[-1]))
        } else if (fx[length(fx)]==0) { # make right grid endpoint bigger
            FindZero(f, c(x[-length(x)], x[length(x)]+gs/2))
        } else if (min(fx)==1) {         # make grid bigger and finer
            FindZero(f, c(x[1]-gs/2, FinerGrid(x), x[length(x)]+gs/2))
        } else {                        # test rejects for small and large x
            a <- which(fx==0)
            return(x[c(a[1]-1, a[1], a[length(a)], a[length(a)]+1)])
        }
    }

    x <- FindZero(f, x)

    ## sanity check
    if (length(x)>1) stopifnot(sapply(x,f)==c(1,0,0,1))

    Disect <- function(x, f, tol) {
        ## given vector x such that f(x[1]) != f(x[2]), find exact location
        ## where f jumps
        stopifnot(length(x)==2)

        if (max(x)-min(x)<tol) {
            return(x)
        } else {
            x <- c(x[1], mean(x), x[2])
            if (f(x[2])==f(x[1])) {    # jump is between x[2] and x[3]
                Disect(x[-1], f, tol)
            } else {                    # jump is between x[1] and x[2]
                Disect(x[1:2], f, tol)
            }
        }
    }

    if (length(x) == 1) {               # this happens if x=maxlength
        return(c(0,x))
    } else {
        return(c(mean(Disect(x[1:2],f,tol)), mean(Disect(x[3:4],f,tol))))
    }
}
